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Abstract. Organisms modulate their fitness in heterogeneous environments by dispers- 

_ ing. Prior work shows that there is selection against "unconditional" dispersal in spatially 

J^ heterogeneous environments. "Unconditional" means individuals disperse at a rate indepen- 

^^ dent of their location. We prove that if within-patch fitness varies spatially and between 

3 two values temporally, then there is selection for unconditional dispersal: any evolutionarily 

^V stable strategy (ESS) or evolutionarily stable coalition (ESC) includes a dispersive pheno- 

(N 

type. Moreover, at this ESS or ESC, there is at least one sink patch (i.e. geometric mean 

pL^ of fitness less than one) and no sources patches (i.e. geometric mean of fitness greater than 

_T one). These results coupled with simulations suggest that spatial-temporal heterogeneity 

• i-H 

r^ due to abiotic forcing result in either an ESS with a dispersive phenotype or an ESC with 

O^ sedentary and dispersive phenotypes. In contrast, spatial-temporal heterogeneity due to 

, biotic interactions can select for higher dispersal rates that ultimately spatially synchronize 

population dynamics. 



> 

(N 

1. Introduction 

O 

O 1 Plants and animals often live in landscapes where environmental conditions vary in space 



> 



2 and time. These environmental conditions may include abiotic factors such as light, space, 



k> 3 and nutrient availability or biotic factors such as prey, competitors, and predators. Since 

c^ 4 the fecundity and survivorship of an individual depends on these factors, an organism may 

5 decrease or increase its fitness by dispersing across the environment. Understanding how 

6 natural selection acts on dispersal in heterogeneous environments has been the focus of 

7 much theoretical and empirical work [lH2T]. 

8 In spatially heterogeneous environments, the general consensus is that there is selection 

9 against unconditional dispersal [31 [H [121 |T6l [20] . Here, unconditional refers to the assumption 

10 that individuals disperse at a rate independent of their location. Using reaction-diffusion 

11 equations, Dockery et al. [12] proved that for two competing populations only differing 

1 
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12 in their diffusion constant, the population with the larger diffusion constant is excluded. 

13 In [ini EU], similar results were proven for populations with non-overlapping generations 

14 living in a patchy environment. Alternatively, in temporally but not spatially heterogenous 

15 environments, Hutson et al. [15J proved that dispersal rates are a selectively neutral trait 

16 for reaction-diffusion models and, thereby, confirmed the numerical observations of McPeek 

17 and Holt [8J for discrete-time, two-patch models. These results imply that the slightest cost 

18 of dispersal would result in selection against dispersal in purely temporally heterogenous 

19 environments. 

20 When there is a mixture of spatial and temporal heterogeneity, the interaction between 

21 competing dispersal phenotypes becomes more subtle. Numerically simulating discrete-time, 

22 two-patch models, McPeek and Holt [8], Parvinen [20], and Mathias et al. [13] found that 

23 more dispersive phenotypes could displace more sedentary phenotypes for certain forms 

24 of spatial-temporal heterogeneity, while evolutionarily stable coalitions of sedentary and 

25 dispersal phenotypes are possible for other forms of spatial-temporal heterogeneity. Hutson et 

26 al. [15j proved that similar phenomena could occur for reaction diffusion equations. However, 

27 analytical criteria distinguishing these outcomes remain elusive. 

28 In this article, we consider the evolution of dispersal for a general class of multi-patch 

29 difference equations varying periodically in time. This periodic variation can be either due 

30 to biotic interactions or abiotic forcing. Our main goals are to analytically identify potential 

31 evolutionarily stable strategies or coalitions for dispersal, characterize the spatial-temporal 

32 patterns of fitness generated by populations playing these dispersal strategies, and use our 

33 results to compare evolutionary outcomes for oscillations due to abiotic forcing versus oscil- 

34 lations due to biotic interactions. 

35 2. Models and Assumptions 

36 To understand the formation of evolutionarily stable coalitions of dispersive phenotypes, 

37 we consider a population consisting of m phenotypes dispersing in an environment consist- 

38 ing of n patches. Let xjit) denote the abundance of phenotype i in patch j. The fitness 

39 of an individual in patch j is assumed to be of the form f-^(t,J2^i^l(t))- ^^ particular, 

40 this assumption implies that phenotypes only differ demographically in their propensity to 

41 disperse. Moreover, we assume that f^{t, ■) is of period p. This periodicity may arise from 

42 exogenous forcing or biological interactions (e.g. over compensating density dependence or 

43 predator-prey interactions). While we do not explicitly model interactions with other species. 
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44 our formulation is sufficiently general to be viewed as the dynamics of a particular species 

45 embedded within a web of interacting species. 

46 We assume that the fraction of phenotype i dispersing from any given patch is di. Of 

47 the individuals dispersing from patch j, a fraction S^j go to patch k. We call S = (Skj) 

48 the dispersal matrix and it characterizes how dispersing individuals are redistributed across 

49 the landscape. Under these assumptions, the interacting phenotypes exhibit the following 

50 population dynamics: 

1=1 J k=l \ 1=1 J 

51 To express this model more succinctly, let x-' = {x{, . . . ,a;^) be the vector of abundances 

52 of the m phenotypes in patch j, Xj = {x} , . . . , x^)"^ (where '^ denotes transpose) be the 

53 vector of abundances of strategy i across the n patches, and Hx-'H = Yl^i^i denote the 

54 total abundance of individuals in patch j. Let x be the matrix with entries x^, F(t, x) be 

55 a diagonal matrix whose j-th diagonal element is /•'(t, Hx-^H), and S{di) = (1 — di)I + diS 

56 where I is the identity matrix. With this notation, the model is represented more succinctly 

57 as 

(2.1) Xi(t + l) = S(c?i)F(t,x(t))x,(t) 2 = l,...,m. 

58 About this model, we make three assumptions throughout this manuscript. 

59 Al: The dispersal matrix S is irreducible and column stochastic (i.e. S has nonnegative 

60 entries, and the entries of each column sum up to one). 



61 A2: (2.1) has a positive period-p point x(l), . . . ,x(p) i.e. ||xi(t)|| > and ||x-'(t)|| > 



62 for all i, j, t, and 

p 
(2.2) J]S(rfi)F(t,x(t))xi(l) = Xi(l) for all 



t=i 



63 A3: F(t, x) is continuous in x. 

64 Assumption Al ensures that individuals or their decedents can move from any patch to any 

65 other patch after sufficiently many generations and there is no direct cost to dispersal. As- 

66 sumption A2 implies that the phenotypes are coexisting in a periodic fashion and occupying 

67 all the patches. Assumption A3 is a basic regularity assumption met by most models. 
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68 3. Main results 

69 We are primarily interested in understanding when a periodically fluctuating collection of 

70 phenotypes can not be invaded by any other phenotype, and what are the spatial-temporal 

71 patterns of fitness at these potential evolutionary end states. To state our main results 

72 precisely, we need two sets of definitions. 

73 Invasion rates, Nash equilibria and evolutionary stability. Let d = {di,. . . ,dm) the 

74 coalition of strategies played by the resident population. If we added a "mutant" phenotype 

75 with dispersal rate d into the population and this mutant population y = (y^, . . . , y^) is very 

76 small, then the resident's population dynamics are initially barely influenced by the mutant's 

77 population dynamics and the mutant's dynamics are approximately given by a linear model 

y(t + l) = S(J)F(t,x(t))y(t). 

78 By standard linearization theorems (see, e.g., p2]), this approximation is valid when the size 

79 of the mutant population is small and the periodic point in A2 is hyperbolic. 

80 The initial fate of the mutant population depends on the invasion rate of strategy d against 

81 the resident population playing strategies d: 

^ i/p 
I{d-d) = g[l[S{d)F{t,^t))^ 

82 where x(t) has period p and ^(A) corresponds to the largest eigenvalue for a non-negative 

83 matrix A. If the invasion rate X(d; d) is greater than one, then the mutant population 

84 grows when its size is small. The ultimate fate of the mutant and resident after the mutant 

85 increases depends on the details of full non-linear model of the resident and invader dynamics. 

86 In particular, following an invasion the asymptotic dynamics in general may no longer be 

87 periodic (i.e. satisfy A2). There are many cases where post-invasion dynamics will remain 

88 periodic (e.g. periodically forced competitive systems as discussed in section 4.1 or when 

89 there is attractor inheritance for a sufficiently small mutation [23]). 

90 If the invasion rate X(d; d) is less than one, then the mutant population declines exponen- 

91 tially when rare and it can not invade. Finally, if X(d; d) = 1, then a mutant may increase 

92 or decrease when rare depending on the details of the nonlinearities of the full model. How- 

93 ever, if it increases, then it does so at a sub exponential rate and, therefore, may be highly 

94 vulnerable to stochastic extinction. An important consequence of our assumption A2 is that 
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95 the invasion rate of mutants with the same strategy as a resident equals one i.e. I(d; d) = 1 

96 whenever d = di for some i. 

97 Using the invasion rates of mutant strategies, we can define several concepts associated 

98 with evolutionary stability. A coalition of strategies d with m > 1 is a mixed Nash equilibrium 

99 provided that 

(3.3) X(d; d) < 1 for all d E [0, 1] 

100 In other words, a mixed Nash equilibrium is a set of strategies in which all mutant strategies 

101 are unlikely to invade due to vulnerability to stochastic extinction. When m = 1, a strategy 



102 satisfying (3.3) is called simply a Nash equilibrium. Under the stronger assumption that rare 

103 mutants decline exponentially to extinction (i.e. X(d; d) < 1 for all d ^ {di, . . . , dm}), d is 

104 an evolutionarily stable coalition (ESC) if m > 1 or a.n evolutionarily stable strategy (ESS) 

105 if m = 1 [21]. More generally, every ESS (respectively, ESC) is a (mixed) Nash equilibirum. 

106 Sources, sinks, and balanced patches. PuUiam [6] introduced the notion of sources and 

107 sink patches for a population at equilibrium. In source patches birth rates exceed death rates, 

108 while in sink patches death rates exceed birth rates. Here, we extend Pulliam's definition to 

109 population exhibiting periodic fluctuations in abundance. A patch is a source if births exceed 

110 deaths "on average" across years, while a patch is a sink if deaths exceed births "on average" 

111 across years. For fitnesses varying in time, the appropriate "average" is the geometric mean: 

^ i/p 

,t=i 

112 If fj < 1, then patch j is a sink. If fj > 1, then patch j is a source. Following McPeek and 

113 Holt [8j, we say patch j is balanced if fj = 1. Individuals remaining in a balanced patch, on 

114 average exactly replace themselves. 

115 Main Results. We have two main results. Our first result implies that if there is spatial 

116 heterogeneity in the within-patch fitnesses under equilibrium conditions, then no coalition of 

117 distinctive pheno types can coexist, there is selection for slower dispersers, and the landscape 

118 supports source and sink patches. In particular, this result implies that for populations 

119 at equilibrium and playing a Nash equilibrium, all patches must be balanced. This result 

120 follows from [16] and generalizes earlier work of Hasings [3] and Parvinen [20j by allowing 

121 for non-diffusive patterns of dispersal (i.e. S need not be symmetric). 
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122 Proposition 1. If p = 1 and F(l,x(l)) is not scalar (i.e. there is spatial heterogeneity), 

123 then all the d.i are equal and positive, X(d; d) > 1 for all < d < di, and at least one patch 

124 is a sink and at least one patch is a source. 

125 Proof. Let A be the diagonal matrix F(l,x(l)) with diagonal entries Ajj = /■'(I, ||x-'(l)||). 

126 Since A is not scalar, Theorem 3.1 in [16] implies that g{S{d)A) is a strictly decreasing 

127 function of d G [0, 1]. Assumption A2 implies that g{S{di)A) = 1 for all i. Hence, all of the 

128 di are equal and X(d; d) > 1 for all < d < di. The di can not equal zero as this would imply 

129 that F(l,Xj(l)) is the identity matrix contradicting the assumption that it is non-scalar. 

130 Since S{di) is column stochastic and A is non-scalar, maxj Ajj > g{S{di)A) = 1 > mini An. 

131 Hence, there is a source patch and a sink patch. 

132 H 

133 Our second result concerns period 2 environments. In contrast to populations at equilib- 

134 rium, this result implies that any ESS or ESC (or more generally, Nash equilibrium) includes 

135 a dispersive phenotype, supports at least one sink patch and supports no source patches. 

136 We are able to prove this result only under the restriction that the dispersal matrix S is 

137 diagonally similar to a symmetric matrix: there exists an invertible diagonal matrix D such 

138 that DSD~^ is a symmetric matrix. This allows for a diversity of movement patterns in- 

139 eluding diffusive movement (i.e. symmetric S) and any form of local movement along a 

140 one-dimensional gradient (i.e. tridiagonal S). A proof is given in Appendix A. It is an open 

141 problem whether this result extends to all irreducible stochastic matrices S. 

142 Theorem 2. Ifp = 2, S is diagonally similar to a symmetric matrix, F(2, x(2)) ^ F(l, x(l)) 

143 (i.e. there is temporal heterogeneity), and F(t,x(t)) is non- scalar for some t (i.e. there is 

144 some spatial heterogeneity), then 

145 (i) // maxj di = 0, then X(d; d) > 1 for all < d < 1. 

146 (ii) If d is a (possibly mixed) Nash equilibrium, then maxj d^ > 0, the set Sinks C 

147 {1,. . . ,n} of sink patches is non-empty, and the remaining patches {1, . . . ,n}\Sinks 

148 are balanced. 

149 4. Applications 

150 To illustrate the utility of our results, we present two applications. The first application 

151 considers populations with compensating density dependence in a periodically forced envi- 

152 ronment. We use Theorem [2] and its proof to determine under what conditions there is an 
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153 ESC of sedentary and dispersive phenotypes. The second application considers populations 

154 with sufficiently overcompensatory density dependence to create oscillatory dynamics. We 

155 use the proof of Theorem [2] to illustrate how the evolution of higher dispersal rates can 

156 synchronize initially asynchronous population dynamics. 

157 4.1. Evolution of dispersal dimorphisms in periodic environments. We consider 

158 competing dispersive phenotypes whose within patch dynamics are given by a periodically 

159 forced Beverton-Holt model. For simplicity, we assume that only the intrinsic fitness A^ of 

160 an individual living in patch j varies in time and space: 



(4.4) f(t,||x^ 



H 



l + a||xJ|| 

161 where a measures the intensity of competition within a patch. Also for simplicity, we assume 

162 that Skj = l/n for all j, k. In other words, dispersing individuals are uniformly distributed 

163 across the landscape. 

164 When there is only one dispersive phenotype (i.e. -m = 1), equation (2.1[) is a sublinear 



165 monotone map (see, e.g., [M] for definitions). Consequently, [25| Thm. 6.1] implies that the 

166 fate of the population depends on the linearization of the system at the extinction state. 

167 The dominant eigenvalue associated with this linearization is given by 

^ i/p 

A = ^( J]S(rfi)F(t,0)' 

a=i 

168 If A < 1, then population goes deterministically toward extinction i.e. x(t) converges to for 

169 all x(0) > 0. Alternatively, if A > 1, then the populations increases when rare and ultimately 

170 converges to a periodic orbit. More precisely, there exists a periodic orbit, {x(l), . . . ,x(p)} 

171 with x(t) ^ for all t, such that x(t) converges to this periodic orbit whenever x(0) ^ 0. 

172 A sufficient condition ensuring A > 1 is 

(4.5) n ^0 ^ ^ ^°^ ^^^ ^' 

173 In other words, all the patches can support a population in the absence of immigration. 



174 While this condition is stronger than what is necessary, for simplicity, we assume that (4.5) 

175 holds for the remainder of this section. 

176 For period two environments where the intrinsic ffiness vary in space and time. Theorem |2] 

177 implies that a sedentary strategy (i.e. di = 0) is not a Nash equilibrium as it can be invaded 
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Figure 1. Nash equilibria for a Beverton-Holt two patch model with A^ = 
2 + ei and A^ = 2 — e^. The strategy (ii = 1 is a local ESS in the regions denoted 
"fast". The mixed strategy d = (1,0) is an ESC in the regions denoted 
"dimorphism". For these regions, the sedentary population only resides in 
the indicated patch. Shading corresponds to the average abundance of the 
sedentary strategy with lighter colors corresponding to higher abundance. At 
the dashed line, all (coalitions of) strategies are a (mixed) Nash equilibrium. 
*s refer to parameter values for which pairwise invasibility plots are shown in 
Fig.|2] 



178 by more dispersive phenotypes. In contrast, for a fully dispersive phenotype (i.e. di = 1), 

179 there is a periodic point x(t) given by 



xiit) 



A2A1 - 1 
a(l + Ail 
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180 where Aj = ^ Yll=i K is the spatial average of the intrinsic fitnesses. Along this periodic 

181 orbit, a computation reveals that the within-patch fitnesses satisfy 

2 2 

nf(t,iix^(t)ii)=ni*. 
t=i t=i * 

182 Thus, Theorem [2] implies that a necessary condition for (ii = 1 to be a Nash equilibrium is 

2 2 

(4.6) l[Xl<l[XtioTallj 

t=i t=i 

183 with a strict inequality for at least one j. This condition for a Nash equilibrium requires 

184 that the geometric mean of the fitness within each patch is no greater than geometric mean 

185 of the spatially averaged fitness. 



186 To illustrate the utility of (4.6), consider an environment where the fitness in each patch 

187 fluctuates between a low value Xbad in "bad" years and a higher value \good in "good" years 

188 i.e. A^ G {\good, Xhad} , Xbad < Xgood, and A^ 7^ A^,^^. To ensure that all dispersal phenotypes 

189 can persist, we assume that the geometric mean ^/Xgood^bad is greater than one. Provided 

190 that there is some spatial asynchrony (i.e. A^ 7^ A^ for some k,j), the necessary condition 



191 (4.6) for a Nash equilibrium of highly dispersive phenotypes simplifies to 

y XgoodXbad ^ -^[Xgood + Afcarf). 

192 This inequality holds strictly as the geometric mean is less than the arithmetic mean. Fur- 

193 thermore, a computation reveals that the geometric mean of fitness within patch j satisfies 



194 for all patches j. Hence, Theorem p^ in Appendix A implies that X(l; rf) < 1 for all (i G [0, 1). 

195 Therefore, for this environment, a highly dispersive phenotype [di = 1) always is an ESS 

196 and all patches are sinks for populations playing this ESS. 



197 When the necessary condition (4.6) for a highly dispersive phenotype to be a Nash equilib- 



198 rium is not met, sedentary dispersers can invade the patches whose average fitness \/nt=i A* 



199 exceeds the spatially averaged fitness y nt=i At- To understand the implications of this in- 

200 vasion, we consider a two patch environment where environmental fluctuations are spatially 

201 asynchronous. More precisely, we assume that spatial average of fitness does not vary in 

202 time (i.e. Xt = X) and the temporal fluctuations et G (—A, A) around this spatial average are 
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Figure 2. Pairwise invasibility plots (PIPs) for a two patch Beverton-Holt 
model with }^^ = 2 + e^ and A^ = 2 — ej. The horizontal and vertical axes 
correspond to resident d\ and mutant d dispersal rates, respectively. The dark 
lines, shaded regions, and unshaded regions correspond to where X(rfi, d) = 1, 

I(cii, (i) < 1, and X((ii, c?) > 1, respectively. In (a) and (b), (ii = is an evolu- 
tionary branching point. In (b), di ~ 0.5 is a convergently unstable singular 
point. In (c), c?i = 1 is convergently stable and evolutionarily stable. In (d), 
di = 1 is convergently stable and a local ESS, but invasable by sufficiently 
sedentary phenotypes. 



203 asynchronous i.e. A^ = A + e^ and A^ = A — e^ for t = 1, 2. The necessary condition (4.6) for 



204 a highly dispersive phenotype to be a Nash equilibrium becomes 



(4.7) 



Alei + €2! < |eie2| and eie2 < 0. 
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205 When (4.7) is not meet, extensive numerical simulations suggest that after the successful in- 

206 vasion of the sedentary dispersers, the populations approach a period 2 orbit x(t). Moreover, 

207 these simulations suggest that d = (1, 0) is an ESC and, consequently, a potential evolution- 

208 ary end point consisting of a dimorphism of sedentary and highly dispersive phenotypes. At 

209 this ESC, one patch is balanced and occupied by both phenotypes, while the other patch is 



210 a sink and only occupied by the dispersive phenotype. (4.7) implies that the ESC occurs 

211 when the temporal correlations of within patch fitness are not too negative (i.e. ei is not 

212 too close to —62 in Fig. fl]). Pairwise invasibility plots (see, e.g., [20] for a discussion of the 

213 interpretation of these plots and the associated terminology) suggest that these ESCs can 

214 be reached by small mutational steps when d = is a convergently stable branching point 



215 (Fig. [2^,b). On the other hand, when (4.7) is satisfied, the highly dispersive phenotype 



216 may or may not be an ESS in the strict sense (Fig. pp,d). Although numerical simulations 

217 confirm that the highly dispersive phenotyperesists invasion attempts by nearby phenotypes 

218 (i.e. X{l,d) < 1 whenever d is sufficiently close to 1), relatively sedentary phenotypes still 

219 may be able invade (Fig. [2]; 



220 4.2. Evolution of synchronicity. Biotic interactions can generate oscillatory dynamics 

221 and, thereby, temporal variation in fitness. To illustrate the feedbacks between evolution 

222 of dispersal and biotically generated oscillations, we consider an extension of the coupled 

223 Logistic map introduced by Hastings [27j. In this model, the local dynamics are determined 

224 by the Logistic fitness function rx{l — x). A fraction d of all individuals disperses randomly 

225 to all patches i.e. a fraction d/k of individuals from patch j arrives in all patches. Under 

226 these assumptions, the dynamics of a single dispersive phenotype is given by 



J ,^ .^ ... ox d 

^t + 1 , , „ , 

n 

fc=i 



d)rxi{l - 4) + -J2rx^{l - x'l) 



227 We note that in the case of n = 2 patches, Hasting's d corresponds to our d/2. 

228 When r > 3 and there are two patches, Hastings [27] has shown that there is an interval 

229 of dispersal rates between and 1 such that there is an out-of-phase stable period two point 

230 (Figs. [3^,c,d). To apply our results, let m = 1, di = d for which there is a stable out- 

231 of-phase periodic point, x(l),x(2) (an explicit formula for this orbit can be found in the 

232 Appendix of [27]), Sjk = 1/2 for I < j,k <2, and p{t, \\x.^\\) = r(l - \\x.^\\). Along this 

233 out-of-phase periodic orbit, j^j F(t, x(t)) is a scalar matrix. Hence, Theorem, [s] implies that 

234 along this out-of-phase periodic orbit T{di, d) > 1 for all di < d < 1. Hence, there would be 
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Figure 3. Evolution of spatial synchronization in a two-patch Logistic model. 
In (a), (c), (d), orbital bifurcation diagrams for phase difference x^ — x^ are 
shown for for the two-patch Logistic equation. All phase differences along 
attractors are plotted in black. The red dashed curve corresponds to the out- 
of-phase, period two orbit which is stable only when it overlaps the black 

regions. In (b), the contours of X{d;d) are plotted along the out-of-phase 
periodic point whenever it is stable. 



235 selection for higher dispersal rates. For this two patch case, this implication of Theorem |3] 

236 also follows from a proposition of Doebeli and Ruxton [101 Pg- 1740] in which they performed 

237 a direct calculation of the eigenvalues. For 3 < r < 3.7, numerical simulations suggest that 
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Figure 4. Evolution of spatial synchronization for a 40 patch Logistic model. 
In (a), the mean phase difference in spatial abundance (i.e. t^ Y1 



402 ^^j^k 



\x-' 



■x''\) 
for a single population exhibiting a two-cycle is plotted as function of its 
dispersal rate. In (b), the strength of selection ^{d,d) for higher dispersal 
rates is plotted as a function of its dispersal rate. Different lines correspond to 
different percentages of initial asynchrony for sedentary phenotype e.g. 50% 
implies that half the patches at rf = where out of sync with the other set of 
patches. 



238 this selection for higher dispersal rates ultimately results in a dispersal rate that spatially 

239 synchronizes the dynamics (Figs, pk-c) at which point all dispersal rates are Nash equilibria. 

240 However, at higher r values such as r = 3.95 (Fig. mi), the destabilization of the out-of-phase 

241 period 2 point (at di ^ 0.38) results in more complex asynchronous dynamics in which case 

242 our results are not applicable and evolution of dispersal may no longer synchronize the 

243 dynamics. 

244 When there are more than two patches and r > 3, out-of-phase 2 cycles can take on a 

245 greater diversity of forms. In particular, one can divide the landscape into two sets of patches 

246 such that patches are synchronous within each set and asynchronous across sets. Due to this 

247 potential spatial asymmetry in these out-of-phase cycles, we have not been able to show 

248 the geometric mean of fitness equals one in all patches. However, numerical simulations for 

249 n = 40 patches and 3 < r < 3.45, suggest that this does occur. In which case. Theorem |3] 

250 implies that there would selection for higher dispersal rates along these asynchronous cycles. 

251 In fact, numerical simulations for 3 < r < 3.45 show that there would be selection for 

252 higher dispersal rates until the dynamics are spatially synchronized (Fig. 111). Moreover, 
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253 these simulations show, quite intuitively, greater initial asynchrony for the dynamics of the 

254 sedentary phenotype result in a stronger selection gradient (Fig.Eb) and require the evolution 

255 of higher dispersal rates to regionally synchronize the dynamics (Fig. |4^ 

256 5. Discussion 

257 We analyzed the evolution of dispersal in spatially and temporally variable environments. 

258 When there is spatial variation in fitness and within patch fitness varies in time between a 

259 lower and higher value, we proved that any evolutionary stable strategy (ESS) or evolution- 

260 arily stable coalition (ESC) includes a dispersive phenotype. In particular, our results imply 

261 a sedentary population can be always be invaded by more dispersive phenotypes. These 

262 results are particularly remarkable in light of earlier analysis showing generically the only 

263 ESS for spatially heterogenous environments without temporal heterogeneity is a sedentary 

264 phenotype [31 [T21 [ISl [SD]- Hence, our results imply this earlier work on discrete-time models 

265 is not generic: arbitrarily small periodic perturbations result in selection for dispersal. 

266 Our results extend a result of Parvinen [20] who proved sedentary phenotypes could be 

267 invaded by phenotypes that disperse uniformly to all patches in every generation. Moreover, 

268 they are consistent with the numerical work of Mathias et al. [13] who found that ESCs 

269 always included highly dispersive phenotypes in discrete-time models with random variation 

270 in the vital rates. In contrast, our results are only partially consistent with the analysis of 

271 reaction-diffusion models by Hutson et al. [15]. Hutson et al. proved that certain forms of 

272 spatio-temporal heterogeneities select for the higher dispersal rate. However, if the frequency 

273 of spatial oscillations (i.e. spatial variation) is too large or too small, then phenotypes with 

274 higher dispersal rates are driven to extinction. A partial explanation for this discrepancy 

275 is the continuity of the per-capita growth rates in the reaction-diffusion models results in 

276 positive correlations in time which may select for slower dispersal rates. 

277 Our analysis and numerical simulations suggest that there are two evolutionary end states 

278 for environments where spatial-temporal heterogeneity is generated by abiotic periodic forc- 

279 ing: either an ESS consisting of a highly dispersive phenotype or an ESC consisting of a 

280 highly dispersive phenotype and a sedentary phenotype. This is partially consistent with the 

281 numerical work of Mathias et al. [13] who always found ESCs of high and low dispersal phe- 

282 notypes. Mathias et al. found that these ESC always could be achieved by small mutational 

283 steps leading to an intermediate phenotype (a convergent singular strategy) at which coali- 

284 tions of faster and slower dispersers could invade and coexist (evolutionary branching). Our 
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285 numerical and analytic results differ from these conclusions in two ways. First, evolutionary 

286 branching occurs at sedentary phenotypes not intermediate phenotypes. Consequently, prior 

287 to the branching event, one of the strategies in the ESC is already present. Second, while 

288 evolution by small mutational states may culminate in a highly dispersive phenotype, an 

289 ESC of highly dispersive and sedentary phenotypes may arise by the invasion of sufficiently 

290 slow dispersers i.e. large mutational steps are required to realize the ESC. 

291 When spatial and temporal heterogeneity is created purely by biotic interactions and 

292 initial conditions (e.g. the coupled Logistic map), we have shown there can be selection 

293 for higher dispersal rates that ultimately may synchronize the population dynamics across 

294 space. Following this synchronization event, dispersal would become selectively neutral. Our 

295 analytic results extend Doebeli and Ruxton's ^0\ two-patch analysis to multiple patches. In 

296 a multispecies context, Dercole et al. [21] found a related result. Numerical simulations of 

297 tritrophic communities with chaotic dynamics showed that evolution of dispersal often drove 

298 these spatial networks to weak forms of synchronization. 

299 Our analysis states that an ESS or a ESC for dispersal result in some patches being sinks 

300 (i.e. geometric mean of fitness less than one) and the remaining patches being balanced 

301 (i.e. geometric mean of fitness equal to one). In the case of ESCs of highly dispersive and 

302 sedentary phenotypes, the sedentary phenotypes only reside in the balanced patches. In the 

303 case of ESS of highly dispersive phenotypes, all patches may be sinks despite the population 

304 persisting. In contrast, evolutionarily stable strategies for conditional dispersal in purely spa- 

305 tially heterogenous environments result in all occupied patches being balanced [T6l ITS] and, 

306 thereby, exhibiting an ideal and free distribution [Ij. By appropriately modifying empirical 

307 methods for distinguishing between source-sink dynamics and balanced-dispersal [HI EH], one 

308 might be able to find empirical support for these alternative, evolutionarily stable, spatial- 

309 temporal patterns of fitness. 

310 Despite extensive progress, there are many mathematical challenges to overcome if we want 

311 to have a better analytical understanding of the evolution of dispersal. Two challenges, of 

312 particular interest here, are going beyond the assumption period 2 environments, and ac- 

313 counting for various forms of costs for dispersal. While period 2 environments can be viewed 

314 as a crude cartoon of seasonal environments, most environmental fluctuations exhibit multi- 

315 pie modes in their Fourier decomposition and have a significant stochastic component [30] . 

316 Despite some progress [ST] |32] , a detailed analytic understanding of the interactive effects 

317 of this environmental stochasticity and dispersal on regional fitness remains largely elusive. 
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318 Alternatively, the ability to disperse or the act of dispersing often is accompanied by costs 

319 to the individual. Dispersing individuals may die before reaching their destination. Alter- 

320 natively, there may trade-offs between dispersal ability and competitive abilities [55| 151]. If 

321 these costs or tradeoffs are sufficiently strong, they can substantially alter the predictions 

322 presented here. For example, Parvinen [20] has shown, quite intuitively, that if costs to dis- 

323 persal are sufficiently strong, there is always selection against dispersal. However, we need 

324 the development of new mathematical methods to unravel the implications of intermediate 

325 costs on the evolution of dispersal in environments with spatio-temporal variation. 

326 Appendix A. Proof of Theorem [2] 

327 The proof of Theorem [2] depends on the following key result which is proven in Appendix 

328 B. We do not impose the assumption that S is irreducible in this result. 

Theorem 3. Let D = diag (di, . . . , (i„) with di,...,dn G (0, oo), and let S be an n x n 
column stochastic matrix such that RSR~^ is symmetric for some diagonal matrix R. For 
t G [0, 1], denote by g{F{t)) the Perron (largest) eigenvalue of 

F{t) = D[{1 - t)In + tS]D-'\{l - t)In + tS]. 

329 Then g{F{t)) is either an increasing function on [0,1] or a constant function on [0,1]. The 

330 latter case happens if and only if D and S commute, equivalently, there is a permutation 

331 matrix P such that PSP^ = Si ® S2 ® ■ ■ ■ ® Sk and PDP^ = dJn^ © d2ln2 ® ■■■ ® 44fe; 

332 where Sj G M„ , for j = 1, . . . ,k, and ni + ■ ■ ■ + Uk = n. 

333 Remark. The above result covers the case of symmetric S. It also covers the case when 

334 S is an irreducible tridiagonal column stochastic matrix; one can use a simple continuity 

335 argument to extend the result to reducible tridiagonal column stochastic matrices. 

336 To prove the first assertion of Theorem |2| assume that maxj di = 0. Assumption A2 

337 implies that F(2,x(2))F(l,x(l)) = /. Theorem js] implies X(d; d) > 1 for all d G (0, 1]. 

338 To prove the second assertion of Theorem pi assume that d is a (possibly mixed) Nash 

339 equilibrium. The first assertion of the Theorem implies that maxj di > 0. Next, suppose 

340 to the contrary that there exists j such that /■'(I, ||x-'(l)||)/-'(2, ||x-'(2)||) > 1 i.e. there is a 

341 source patch. Then 

X(d;0) = maxv/p(l,||x^(l)||)p-(2,||x^X2)||) > 1 
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342 and by continuity X(d; d) > 1 for all d > sufficiently small. As this contradicts the 

343 assumption that d is a Nash coalition, it follows that /■'(I, ||x''(l)||)/-'(2, ||x-'(2)||) < 1 for all 

344 j. Finally, suppose to the contrary that /•'(I, ||x''(l)||)/''(2, ||x-'(2)||) = 1 for all j. Theorem|3] 

345 implies that X(d; maxj di) > 1 which contradicts the fact that X(d; maxj di) = 1. 

346 Appendix B. Proof of Theorem [3] 

347 Denote by \\A\\ the operator norm of the matrix A. The proof of Theorem Is] depends on 

348 the following. 

349 Theorem 4. Suppose A G M„ is nonzero and satisfies || J+y4|| > 1. Then || J+M|| > ||/+y4|| 

350 for all t > 1 and one of the following condition holds. 

351 (a) The function t h-> ||/ + tA|| is increasing for t > 1. 

352 (b) There is a unitary matrix U such that U*AU = 0^ © A, where A G M^-k is invertible 

353 and satisfies \\I„-k + ^|| < 1- Consequently, there is t* > 1 such that ||/„,_fc + t*A\\ = 1 and 

354 the function t \-^ \\I + tA\\ has constant value 1 for t G [1, t*] and is increasing for t > t*. 

Proof. Let -u be a unit vector such that ||/ + A\\ = ||(/ + A)u\\ and Au = au + /3v, where 
{u, v} is an orthonormal family. By our assumption, 

\\(I + A)u\\ = \\{1 + a)u + f3v\\ = |l + «p + |/3|2>l, 

i.e., 

2Re(a) + lap + |/3|2 > 0. 

355 Thus, for t > 1, 

||/ + tA|| > \\{I + tA)u\\ 

= |l + to|^ + |t/3|^ 

= |1 + a|2 + |/3|2 + 2(t - l)Re(a) + (t^ - l)[\a\^ + |/3|2] 

= ||/ + A\\ + (t - l)[2Re(a) + \a\^ + \/3\^] + tit - l)[\a\^ + \/3\^] 

> \\i + M- 

356 (a) Suppose there is unit vector u in the above calculation satisfying Au ^ 0, i.e., (a, /3) ^ 

357 (0,0). Then the last inequality in the calculation is a strict inequality. Thus, ||/ + tv4|| > 

358 ||/ + A||. 
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Now, for any t^ > 1, we have ||/ + tQA\\ > \\I + A\\ > 1. If mq is a unit vector satisfying 
II (/ + to^)'Wo|| = 11-^ + ^0^11 5 then tQAuo ^ 0; otherwise, ||/ + to^ll = 1- Consequently, if we 
replace A by t^A in the above proof, we have 

||/ + t(M)|| > II/ + MII 

359 for any t > 1. Thus the function t i— t- ||/ + ty4|| is increasing for t > 1. 

(b) Suppose Au = for every unit vector u satisfying ||(/ + A)-u|| = ||/ + A||. In particular, 
we have 

||/ + A|| = \\{I + A)u\\ = \\u\\ = 1. 

Let U be unitary such that U*AU is lower triangular form with the first k diagonal entries 
equal to zero, and the last n — k diagonal entries nonzero. Since 

||/ + A|| = \\I + U*AU\\ = 1, 

we see that 

e*(/ + U*AUy{I + U*AU)ej = ||(/ + U*AU)ejf < 1 

360 for j = 1, . . . , k. (Here {ei, . . . , e„} is the standard basis for (D".) As a result, we see that 

361 U*AU = Ofc © A, where A e M„_fe is invertible. 

362 Since Au = for every unit vector u satisfying ||(/ + A)u\\ = \\I + yl|| = 1, we conclude 

363 that ||(/„_fc + A)v\\ < 1 for any unit vector v E (D""^. Thus, ||/„_fc + A\\ < 1. 

364 Note that A ^ implies A is non-trivial, i.e., k ^ n. For sufficiently large t, we have 

365 ||/n-A; + tA\\ > 1. Let t* be the smallest real number in (1, cxd) such that ||/n-fc + t*A\\ = 1. 

366 Since A is invertible, case (a) must hold and the function t t— ;■ ||/n-i + tA\\ is increasing for 

367 t > t*. Hence for t > t*, we have ||/ + tA\\ = \\In-k + A\\ and the function t \-^ \\I + tA\\ is 

368 increasing. ■ 

369 Proof of Theorem^ Note that F{t) = D[{1 - t)/„ + tS]D-'^[{l - t)/„ + tS] and 

F{t) = RF{t)R-^ 

= RDR-^ 1(1 - t)In + tRSR-^]RD~^R~^ [(1 - t)/„ + tRSR-^] 
= D[{1 - t)In + tRSR-^]D-\l - t)In + tRSR-^] 

have the same eigenvalues and hence the same spectral radius. Here, we use the fact that 
RDR^^ = D a.s Ris a diagonal matrix. So, S = RSR~^ is symmetric with largest eigenvalue 
equal to 1, and all other eigenvalues lying in [—1, 1]. Moreover, if 
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370 then D-^/^F{t)D^'^ = B{t)B{tf so that ^(F(t)) = g {D-^/^F{t)D^'^) = ||5(t)f . 

Suppose < to < to + t < 1. Let Ao = tQD^''^{S - In)D-^l'^. Then 

\\I + A4 = \\B{t,)\\>g{B{to)) = l. 

371 By Theorem |4| 

(B.8) ||5(to + t)|| = ||/+(l + t/to)^o|| > ||/ + Ao|| = ||5(to)||. 

372 Thus, g{F(t)) = ||i?(t)|p is non-decreasing on [0, 1]. Moreover, by TheoreniHthe inequahty 



373 in (B.8 ) is an equahty if and only if Aq is unitarily similar to 0^ © Aq where Aq is invertible. 



374 Thus, 

375 (1) the null space of Aq has dimension k, and 

376 (2) Aqv = if and only if v^Aq = 0. 

377 From (1), we see that the eigenspace of the Perron root of the matrix D^^'^SD^^^'^ has 

378 dimension k. So, the matrix is permutationally similar to a (/c + 1) x [k + 1) upper triangular 

379 block matrix so that each of the first k diagonal block is square irreducible and has eigenvalue 

380 1, and the last block has eigenvalue less than 1. Since S is symmetric and D is diagonal, there 

381 is a permutation matrix P such that PSP^ = Si © ■ ■ ■ © 5*^+1 and PDP^^ = -Di © ■ ■ ■ © -Djt+i 

~ 1/2 1/2 

382 with Sj, Dj G Mn^ and ni + ■ ■ ■ + n^+i = n. But then Dfj_^^Sk+iD^j^^ cannot have spectral 

383 radius less than 1. So, 5^+1 must be vacuous. 

384 Next, we show that each Dj is a scalar matrix. To this end, note that we can construct a 

385 null vector of PAqP^ by extending a Perron vector ui G IR"^ of 5*1 to a vector vi in H" by 

386 adding n — ni zeros. Then PAqP^vi = 0. By (2), we see that v^PAqP'^ = 0. It follows that 

387 D^ SiD-^ ui = ui and u[D^ SiD^ = uj. Equivalently, SiD^ ui = D^ u\ and 

388 uJD-^ Si = ujD^ . Since S^i is symmetric, nonnegative, and irreducible, the eigenspace of 

389 the Perron eigenvalue is one dimensional and there is a positive vector w such that Siw = w 

390 and w'^Si = w'^. Hence D^ mi is a multiple of w and u[D^ is a multiple of w'^ . Hence 

391 Di is a scalar matrix. Similarly, one can prove that D2, . . . , Dk are scalar matrices. Since 

392 PSP'^ = Si © ■ • • © ^fc, we have PSP"^ = Si®---®Sk. 

393 Conversely, if 5* and D commute, then F{t) = {tl + (1 — t)S)(tI + (1 — t)S) is column 

394 stochastic for all t G [0, 1]. So, Q{F{t)) = 1 for all t G [0, 1]. ■ 
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395 Remark Note that the conclusion of Theorem |3] fails if S is not diagonally similar to a 

/O 1 0\ 

396 symmetric matrix. For example, if 5" = I 1 © In-s, and D = diag(3,2, 1) © In-s, 

397 then g{F{0)) = g{F{l)) = 1, but g{F{l/2)) > 1. One may perturb 5* slightly, say, replacing 

398 it by (1 — 6)S + eJ for a small e > so that g{F{t)) is not monotone, where J is the matrix 

399 with all entries equal to 1/n. 
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